Orthogonal measurement-assisted quantum control 
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Existing algorithms for the optimal control of quantum observables are based on locally optimal 
steps in the space of control fields, or as in the case of genetic algorithms, operate on the basis of 
heuristics that do not explicitly take into account details pertaining to the geometry of the search 
space. We present globally efficient algorithms for quantum observable control that follow direct or 
close-to-direct paths in the domain of unitary dynamical propagators, based on partial reconstruction 
of these propagators at successive points along the search trajectory through orthogonal observable 
measurements. These algorithms can be implemented experimentally and offer an alternative to 
the adaptive learning control approach to optimal control experiments (OCE). Their performance 
is compared to that of local gradient-based control optimization. 

PACS numbers: 03.67.-a,02.30.Yy 



I. INTRODUCTION 

The optimal control of quantum dynamics is receiv- 
ing increasing interest due to widespread success in lab- 
oratory and computational experiments across a broad 
scope of systems. With these promising results it be- 
comes imperative to understand the reasons for success 
and to develop more efficient algorithms that can in- 
crease objective yields to higher quality. In the compu- 
tational setting, the expense of iteratively solving the 
Schrodingcr equation necessitates faster algorithms for 
the search over control field space if these methods are 
to be routinely employed for large-scale applications. In 
the laboratory setting, although closed-loop methodolo- 
gies have encountered remarkable success, the search al- 
gorithms currently used do not typically attain yields as 
high as those that can be achieved using computational 
algorithms. 

Recently, significant strides have been made towards 
establishing a foundation for the systematic develop- 
ment of efficient OCT algorithms based on the observa- 
tion that the landscape traversed by search algorithms in 
the optimization of quantum controls is not arbitrarily 
complicated, but rather possesses an analytical struc- 
ture originating in the geometry of quantum mechanics 
[l[ . This structure should allow not only rationalization 
of the comparative successes of previous quantum con- 
trol experiments, but also analytical assessment of the 
comparative efficiencies of new algorithms. 

Prior work established important features of these 
landscapes, in particular, their critical topologies @,||. 
The critical points of a control landscape correspond to 
locally optimal solutions to the control problem. The 
most common objective in quantum optimal control is 
maximization of the expectation value of an observable. 
The landscape corresponding to this problem was shown 
to be almost entirely devoid of local traps, i.e., the vast 
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majority of local suboptima are saddles, facilitating the 
convergence of local search algorithms. The number of 
local suboptima, as well as the volumes of these criti- 
cal regions were calculated. However, the relationship 
between the topology of quantum control landscapes 
and their geometry, which would dictate the behavior 
of global search algorithms, was not explored. 

Thus far, optimal control algorithms for quantum ob- 
servables have not exploited the geometry of quantum 
control landscapes, which may simplify control optimiza- 
tion compared to that for classical systems. Indeed, the 
majority of quantum optimal control algorithms to date 
have aimed at optimizing an objective functional, such 
as the expectation value of an observable operator, di- 
rectly on the domain of time-dependent control fields 
s(t). Typical approaches to quantum control optimiza- 
tion use the information in the measurement of a single 
quantum observable to guide the search for optimal con- 
trols; the simplest approach is to randomly sample single 
observable expectation values at various points over the 
landscape and use genetic algorithms (GA) to update 
the control field. A recent experimental study Q demon- 
strated at least a two-fold improvement in optimization 
efficiency through the use of local gradient algorithms 
rather than GA, but the important question remains as 
to whether global algorithms for quantum control that 
are not "blind" like GA can be implemented in an ex- 
perimental setting. 

When control optimization seeks to optimize the ex- 
pectation value of an observable by following, e.g., local 
gradient information on the domain of controls, the ge- 
ometry of the underlying space of quantum dynamical 
propagators U (N) is not explicitly exploited. In particu- 
lar, optimal control algorithms that are based on locally 
minimizing an objective function on the domain of con- 
trol fields do not follow globally optimal paths in U(N). 
This approach tends to convolute the properties of the 
map between control fields and associated unitary prop- 
agators with the properties of the map between unitary 
propagators and associated values of the objective func- 
tion. 
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An alternative approach to observable maximization 
is to first solve numerically for the set of unitary matrices 
U that maximize the expectation value Tr([/p(O)[/^0) 
of the observable 8, and then to determine a control 
field e(t) that produces that U at time t = T @, ||. 
For restricted Hamiltonians in low dimensions, analyti- 
cal solutions for the optimal control field e(t) have been 
shown to exist. However, to date, numerical algorithms 
for the optimization of unitary propagators in higher 
dimensions have operated solely on the basis of local 
gradient information, such that the global geometry of 
U(N) is again not exploited. 

The variational problems of optimal control theory 
admit two types of minimizers. Denoting the cost func- 
tional by J, according to the chain rule, 

5 J _ dJ SU 
Se(t) ~ dU ' Se(t) ' 

The first type of minimizer corresponds to those control 
Hamiltonians that are critical points of the control ob- 
jective functional, but are not critical points of the map 
between control fields and associated dynamical propa- 
gators (i.e., points at which = 0, while the Frechet 
derivative mapping from the control variation de(t) to 
5U(T) at t = T is surjective). The second type corre- 
sponds to critical points of the latter map (i.e., points 
at which the mapping from Se(t) to SU(T) is not lo- 
cally surjective) |7| 1 . Critical points of the first type, 
which are referred to as kinematic critical points or nor- 
mal extremal controls, are either global optima or sad- 
dle points, but never local traps [2, HJ- Recent work in 
quantum optimal control theory suggests that the crit- 
ical points of the map e(t) — > U(T), called abnormal 
extremal controls, are particularly rare (i.e., there are 
generally fewer critical points compared to classical con- 
trol problems) [l], @] • 

Because the objective function J is a complete func- 
tion of U(T), the maximal achievable optimization ef- 
ficiency is ultimately determined by the properties of 
the map between control fields and unitary propagators, 
e(t) — » U(T). Irrespective of the corresponding observ- 
able expectation value, updating the control field to pro- 
duce a unitary propagator that is close to the current 
propagator will typically be computationally inexpen- 
sive. Therefore, following a direct route in the space 
of unitary propagators is expected to be more efficient 
in quantum control than following a gradient flow on 
the space of objective function values that maps to a 
longer path in U(N). As will be shown, the scarcity 
of critical points of the map e(t) — > U(T) in quantum 
control problems implies that it is surprisingly simple to 
track arbitrary paths in IA{N) during optimization, at 
least for certain families of Hamiltonians. However, it 
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is not uncommon to encounter regions of U{N) where 
numerically, the relevant differential equations are ill- 
conditioned. The ability to selectively avoid such singu- 
lar regions, which correspond to abnormal extremals, is 
desirable. One way to achieve this goal is to constrain 
the search trajectory to only roughly follow a predeter- 
mined path in IA(N). 

In this paper, we develop globally efficient algorithms 
for the optimization of quantum observables that ex- 
ploit the geometry oilA{N) by approximately following a 
predetermined path in the space of quantum dynamical 
propagators. This approach to globally efficient quan- 
tum control optimization is based on making a partial 
tomographic set of measurements at various steps along 
the search trajectory. A complete tomographic set of 
observations is a set that is adequate for the estima- 
tion of all the iV 2 parameters of the unitary propagator 
U(T) @, Q; a partial tomographic set reconstructs only 
a subset of these parameters. The goal of this approach 
is to reap the benefits of unitary matrix tracking with- 
out encountering the associated singularities. As such, 
the approach attempts to leverage the methodologies of 
quantum statistical inference [1 01 ] in order to reduce the 
search effort involved in solving quantum control prob- 
lems. 

These experimentally-implementable algorithms for 
quantum control optimization can be simulated by em- 
ploying a generalization of the diffeomorphic homotopy 
tracking methodology D-MORPH (diffeomorphic modu- 
lation under observable response-preserving homotopy) 
flU [l2l [HI]. In contrast to observable-preserving dif- 
feomorphic tracking, the orthogonal observable track- 
ing algorithm developed and applied here identifies 
parametrized paths e(s,i) that follow a given predeter- 
mined trajectory through U{N). In both denu- 
merably infinite number of solutions exist to the tracking 
differential equations; paths e(s, t) that optimize desir- 
able physical features of the control field can be tuned 
through the choice of an auxiliary free function. We will 
show that a primary difference between scalar and vec- 
tor observable tracking is that the trajectory followed in 
IA(N) in the former case is highly sensitive to changes in 
the system Hamiltonian, whereas the [/-trajectory fol- 
lowed in the latter case can be rendered largely system 
independent by employing a larger set of orthogonal ob- 
servables. This suggests that, besides its usefulness as 
an optimization algorithm, orthogonal observable track- 
ing can reveal universal features underlying the com- 
putational effort involved in quantum optimal control 
searches across diverse systems. 

In addition, we compare the trajectories in U{N) fol- 
lowed by standard OCT gradient-following algorithms 
with those that track optimal paths in the dynamical 
group, for various Hamiltonians, in order to determine 
how the geometry of the underlying space affects the 
convergence of experimental and computational control 
optimizations that exploit only local gradient informa- 
tion. In so doing, we will show that there exists a spe- 
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cial relationship between the gradient flow on e(t) and a 
particular (global) path U (s) in the domain of unitary 
propagators, namely the gradient flow of the objective 
on U(N), which offers insight into the convergence prop- 
erties of the former. 



II. QUANTUM OPTIMAL CONTROL 
GRADIENT FLOWS 



Local algorithms for quantum optimal control, 
whether numerical (OCT) or experimental (OCE), are 
typically based on the gradient of the objective function. 
In this section, we review the properties of the gradient 
for quantum observable expectation value maximization, 
and dissect these properties into system(Hamiltonian)- 
dependent and universal system-independent parts. 

A generic quantum optimal control cost functional can 
be written: 



J = ${U(T),T) 
i T 

Tr 



Re 



A / | e(t) \* dt (1) 
Jo 

where H is the total Hamiltonian, (3(t) is a Lagrange 
multiplier operator constraining the quantum system 
dynamics to obey the Schrodinger equation, e(t) is the 
time-dependent control field, and A weights the impor- 
tance of the penalty on the total field fluence. Solutions 
to the optimal control problem correspond to = 0. 
The functional <!>, which we refer to as the objective 
function, can take various forms. The most common 
form of 4> is the expectation value of an observable of 
the system: 

$(£/) = Tr(C/(T) ( o(0)C/ t (T)e) 

where p(0) is the initial density matrix of the system 
and O is an arbitrary Hermitian observable operator Q . 

An infinitesimal functional change in the Hamiltonian 
5H(t) produces an infinitesimal change in the dynamical 
propagator U(t, 0) as follows: 



SU(t,0) 



U(t,t')6H(t')U(t',0)dt' 



where SH(t) = V e H(t)-Se(t). The corresponding change 
in 4> is then given by 



5$ = - 



ft Jo 



Tr( [9(T), U\t, 0)5H(t)U(t, 0) ]p(0))dt, 



where 6(T) = W(T)p(0)U(T). In the special case of the 
electric dipole approximation, the Hamiltonian assumes 
the form 

H(t) = H - p ■ e(t) 



where Hq is the internal Hamiltonian of the system and 
/j, is its electric dipole operator. fi(t) is given by 

ix{t) ee U\t, 0)fJ,U(t, 0) = ihM(t) 

where M(t) = |t/t (t,0)V E H(t)U(t,0). Within the elec- 
tric dipole approximation, the gradient of <£> is [3|: 



(5$ 



:Tr{[e(T),/x(t)]p(0)} 



Se(t) h 

= l Tt Y J Pm\®(TMt)-n{t)Q{T)\i) 



h ■ 

N 



1=1 3 = 1 



(2) 



where the initial density matrix is given as p(Q) — 

E"=i^N)(*l'Pi > ■■■ > Pn > 0, Z)"=iP< = !• The 
assumption of local surjectivity of e(t) — > U(T) implies 

that the functions (i\n(t)\j) are A^ 2 linearly independent 

functions of time. The functions 



(i\Ui(T)QU(T)\j)(j\U\t)»U(t)\i) - 

(i\tf(t)»U(t)\j){j\uHT)Gtf(T)\i) 



(3) 



therefore constitute natural basis functions for the gra- 
dient on the domain e(t). We are interested in the global 
behavior of the flow trajectories followed by these gradi- 
ents, which are the solutions to the differential equations 



de(s,t) , / w S<^(s,T) 



ds 



Se(s, t) 



(4) 



where s > is a continuous variable parametrizing the 
algorithmic time evolution of the search trajectory, and 
a is an arbitrary positive constant that we will set to 
1. The existence of the natural basis ((3]) indicates that 
these flow trajectories evolve on a low-dimensional sub- 
space of s(t). However, the gradient flow equations 
cannot be integrated analytically for arbitrary internal 
Hamiltonians H , precluding a deeper understanding 
of the global dynamics of the search process. In fact, 
these dynamics do not have universal (Hamiltonian- 
independent) properties. The explicit path followed by 
the search algorithm on s(t) depends on the solution 
to the Schrodinger equation for the particular system 
Hamiltonian and cannot be expressed analytically. 

Because the objective functional $ is explicitly a func- 
tion of U (T) , any universal properties of the global ge- 
ometry of the search dynamics must be investigated on 
this domain. These search dynamics are governed by 
the gradient flow of $ on the domain U(N), given by 

f = V*(*7). 
as 

The tangent space oiU(N) at any element U € U(N) is 
T V U(N) - {I V. <> = -SI, D C NxN }, 
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where f2 is an arbitrary skew-Hermitian matrix, and the 
directional derivative for a function $ defined on U(N) 
is 

D<f>u(UQ) ee Tr ((V$(C/)) t t/0) . 

The directional derivative of the objective functional $ 
along an arbitrary direction UVt in TjjU(N) can then be 
written 

= Tr ([/>((}), f/ + 6Z7]n) 
allowing us to identify the gradient of <E> on U(N) as 

v$ = -?/[p(o), c/tet/j = [e, Up(o)u*]u. 

Therefore, the equations of motion for the gradient flow 
lines of objective functional $ are 



[e, u P (o)u^]u = -up{o)u*eu + eup(o). (5) 



dU 
ds 



In section I VII below, we integrate these equations to 
obtain the trajectories U(s) followed by gradient algo- 
rithms on U(N) over algorithmic time < s < oo. 

The essential question arises as to the relationship 
between the gradient flow on e(t) and that on U{N). 
The gradient on e{t) is related to the gradient on U(N) 
through 



<5$ 



ESUij d<f> 
Se(t) dUa 



(6) 



Now suppose that we have the gradient flow of e(s, t) 
that follows ([4]) and let U (s) be the projected trajectory 
on the unitary group U{N) of system propagators at 
time T, driven by e(s,t). The algorithmic time deriva- 
tive of U (s) is then 



dU tJ (s) 



5Uij(s) de(s, t) 
ds J 5e(s,t) ds 

which, combined with ((3]) and ([5]), gives 



dt 



ds 



T 

Se(s,t) 



P.Q 



5e(s,t) dU pq 



dt. 



(7) 



(8) 



It is convenient to write this equation in vector form, re- 
placing the N x N matrix U (s) with the iV 2 dimensional 
vector u(s): 



du(s) 
ds 



Su(s) Su T (s) 

o Se(s,t) 6e(s,t) 



dt 



V$[u(s)] := 



G[ £ ( S ,t)]V$[u( S )] (9) 



where the superscript T denotes the transpose. Thus 
the projected trajectory from the space of control field 



is different from that driven by the gradient flow in the 
unitary group: 



dU(s) 
ds 



= V$[U(s)]. 



(10) 



This relation implies that the variation of the propaga- 
tor in U(N) caused by tracking the gradient flow in the 
space of control field is Hamiltonian-dependent, where 
the influence of the Hamiltonian is all contained in the 
iV 2 -dimensional symmetric matrix G[e(s,t)]. 



III. UNITARY MATRIX FLOW TRACKING 

The matrix G [e(s, t)} in equation © above indicates 
that the convergence time for local gradient-based OCT 
algorithms may vary greatly as a function of the Hamil- 
tonian of the system. Given the decomposition of the 
gradient into Hamiltonian-dependent and Hamiltonian- 
independent parts, the natural question arises as to 
whether the Hamiltonian-dependent part can be sup- 
pressed to produce an algorithm whose convergence time 
will be (approximately) dictated by that of the unitary 
gradient flow, irrespective of the system Hamiltonian. 

In order for the projected flow from e(t) onto U(T) to 
match the integrated gradient flow on U(T), the quan- 
tity de g S f^ that corresponds to movement in each step 
must satisfy a generalized differential equation: 



dU(s) 



SU(s) de(s,t) 



ds J 5s(s,t) ds 



dt=V<f>[U(s)]. (11) 



In the dipole approximation, this relation becomes the 
following matrix integral equation: 



M (M)^dt = C/t( s )v$[[/( s )] 



where p(s,t) = W (s , t) pXJ (s , t) . When $ is the observ- 
able expectation value objective function, we have 



% (M) ___M dt = _ [ p ,uHs)QU(s)] . 

OS 



On the basis of eigenstates, the matrix integral equation 
is written 

IMi(s, t)^^dt - ih{i\Ul{s, T)V$ [U(s, T)\ \j). 
os 

, , . (12) 

To solve this equation, we first note that the flexibility 
in the choice of the representation of the variation in 
e(s,t) allows us to expand it on the basis of functions 

Pij {s,t) : as 



de(s, t) 
ds 



— y ] XijP'ij ( s j t). 



1,3 
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Inserting this expansion into the above equation pro- 
duces 

y^Xpgis) fj,ij(s,t)fj, pq (s,t)dt = 

ih{i\tf(s,T)v$[U(s,T)}\j). (13) 
If we denote the correlation matrix G(s) as 

Gij.p q (s) = / fXij(s,t)fj. pq (s,t)dt = 
Jo 

f (Ms,t)\j)(p\fi(s,t)\q)dt : (14) 
Jo 

(as in eqn ([9]) above, but now specifically in the case of 
the dipole approximation) and define 

Ai 3 -(s) =ih(i\U*(a,T)s7$[U(s,T)] \j), 

it can be shown (l5| that the matrix integral equation 
(|12|) can be converted into the following nonsingular iV 2 - 
dimcnsional algebraic differential equation: 

g=/ s +( V (A)-a)Vy M (*)) (15) 

where f a = f s (t) is a "free" function resulting from 
the solution of the homogeneous differential equation, 
the operator v vectorizes its matrix argument and a = 

/ T »(/i(t))/.dt. 

Solving this set of A^ 4 scalar differential equations re- 
quires that the iV 2 x N 2 matrix G is invertible. The in- 
vertibility of this matrix is equivalent to the claim that 
the map e(t) — > U(N) between control fields and uni- 
tary propagators is surjective, such that it is possible 
to reach any U(s + 1) infmitcsimally close to U(s) in a 
vanishingly small step. Thus, a necessary condition for 
the existence of a well-determined search direction is the 
full-rank of the Jacobian, i.e., 

ranki^fL = dim[W(JV)] = iV 2 
oe(s,t) 

which is equivalent to the requirement of local surjectiv- 
ity of the map e(t) ->U(T). 

The problem is undetermined because of the rank of 
the matrix is lower than the number of variables to be 
solved, which results from the fact that the optimal con- 
trol problem itself is undetermined with a multiplicity 
of solutions. Each "free function" f s corresponds to a 



unique algorithmic step in e(t); modulating this function 
allows for systematic exploration of the set of functions 
e(s , t) that are compatible with the gradient step on U 

As in the case of the gradient V$[e(t)] , the flow on e(t) 
that tracks the U -gradient can be expressed in terms of 
a set of maximally N 2 linearly independent functions 
of time, but whereas the former is a unique functional 
derivative, the latter is highly degenerate. The former 
are explicitly determined by the functions n(t), while the 
latter are underdetermined by these functions; only N 2 
linearly independent components of d£ Q S ^ are explicitly 
determined by n(t), the rest remaining unspecified. 

Although the gradient step de ^ t ' ) = V$[e(s,i)] is al- 
ways locally the direction of fastest decrease in the ob- 
jective function at e(t), the path e(s,t) derived from 
following this gradient has no universal (Hamiltonian- 
independent) global geometry, since $ is not explicitly 
a function of e(t). It is known [3j that this path will not 
encounter any traps during the search, but beyond this, 
the geometry can be expected to be rugged and globally 
suboptimal. Unlike the gradient V$[e(i)], the algorith- 
mic step above follows the gradient flow on U(N) (in 
the limit of infinitesimally small algorithmic time steps). 
The iV 2 functions /i(s, t) are calculated during the evalu- 
ation of V$[er(i)]; hence, the computational overhead in- 
curred by following this flow corresponds to that needed 
to compute the N A elements of G(s) and invert the ma- 
trix, at each algorithmic time step. This flow respects 
the geometric formulation of the optimal control objec- 
tive function in terms of U(T) rather directly in terms 
of e(t). As we shall show below, the global geometry 
of this path can be completely determined analytically 
for objective function The functions /it(s,£) contain 
all relevant information about the quantum dynamics, 
whereas the functions V$(t/) contain complete informa- 
tion about the geometry of the search space. 

Incidentally, matrix integral equation (I12p can be 
reformulated to provide further insight into the rela- 
tionship between e(i)-gradient and [/(T)-gradient flows. 
Equation (|12p can be rewritten 

(—iH [9(s, T),p], n(s, t)). (16) 

It can be shown that if the time-dependent dipole oper- 
ator fx(s,t) displays the Dirac property 



N N 

i=i j>i 



G 



the corresponding nonsingular initial value problem is 

de(s, t) 



ds 



-ih[&(s,T),p(0)},v(s,t)}, 



which is effectively identical to the e(£)-gradient flow for 
observable functional <£>. As such, the extent to which 
condition (JTSJ) is satisfied for a given Hamiltonian will 
determine the faithfulness with which this flow tracks 
the [/(T)-gradient. 

Of course, a multitude of other flows could be sub- 
stituted for the RHS of equation fT2")) . In section |VH 
we will integrate the [/(T)-gradient flow and show that 
it does not follow a globally optimal path. Since we 
are interested in global optimality, we should choose a 
flow that follows the shortest possible path from the ini- 
tial condition to a unitary matrix that maximizes the 
observable expectation value. It can be shown [j| that 
a continuous manifold of unitary matrices W maximizes 
$(T). These Ws can be determined numerically by stan- 
dard optimization algorithms on the domain of unitary 
propagators [HI]. The shortest length path in U(N) 
between U(0) and an optimal W is then the geodesic 
path that can be parameterized as U(s) = U(0) exp(iAs) 
with A = — i \og(W Jf U (0)) where log denotes the com- 
plex matrix logarithm with eigenvalues chosen to lie on 
the principal branch — n < 8 < n. Thus, if we set 
A I3 (s) = (i\A\j) = (i\ - ilog(W^U(s)\j), the tracking 
algorithm will attempt to follow this geodesic path. Be- 
cause this choice of A does not represent the gradient of 
an objective function, the optimization will not converge 
exponentially to the solution, but rather will continue 
past the target matrix W unless stopped [HI]. On the 
other hand, it is in principle possible to choose A that 
results in the algorithm tracking the same step in U(N), 
but at a rate that depends on algorithmic time s. 2 

Due to the nonlinearity of the differential equations 
above, errors in tracking will inevitably occur, increas- 
ing the length of the search trajectory beyond that of 
the minimal geodesic path (see below). These errors 
will naturally be a function of the system Hamiltonian. 
It is of interest to examine the dependence of matrix 
flow tracking errors on the Hamiltonian by continu- 
ously morphing the Hamiltonian during the optimiza- 
tion. This efficient approach to Hamiltonian sampling 
will allow a more systematic comparison of the efficiency 
of global OCT optimization with that of local gradient- 
based OCT, which is expected to be much more system- 
dependent. 

Hamiltonian morphing can encompass changes in both 
the system's internal Hamiltonian and the dipole op- 
erator. We assume these matrices can be written as 
functions of the algorithmic step s as H(s) — Hq(s) + 
p,(s)e(s,t). Since we have 



2 In the case that the control system evolves on a subgroup of 
U (N), e.g. SU(N), the geodesic on that subgroup can be tracked 
instead. 



dn(s,t) dH Q (s) d/i(s) de(s,t) 

e{s,t) - fi(s)- 



ds ds ds 

we can rewrite eqn (|12j) as 

dU(s) >< 



ds 



dtia {s,t,T) — h 

v OS 

ai (s,t,T)e(s,t) + a 2 (s,t,T)) =0 (17) 



where ao = /i(s,£), a± = dM ^'^ and a-i — d7 ^°„ ■ Thus, 
if we define 



ds 



b(s, T) = — (i/(ai)e(«, t) + u(a 2 ))dt, 
Jo 

where v denotes the iV 2 -dimensional vectorized Hermi- 
tian matrix as above, we can rewrite the matrix integral 
equation (|12[) as 



It 



(s, t)^^dt = - [p, U^s)QU(s) } + b(s, T). 



Therefore in the case of combined Hamiltonian mor- 
phing and unitary tracking, the D-MORPH differential 
equation for the control field becomes 



^ = f s +[v(A) + b(s,T)-a) G-y^*)) (18) 

Even if G is invertible, it is possible that it is nearly 
singular, resulting in large numerical errors during the 
solution to the differential equation. It is convenient 
to assess the nearness to singularity of G by means of 
its condition number C, namely the ratio of its largest 
singular value to its smallest singular value, i.e., if G _1 = 
F[diag(lM)]f/ T ,C=^. 

As mentioned, tracking errors can also originate due 
to the omission of higher order functional derivatives, 

such as $ e (j-}2 , in equation IT7|) . These may in princi- 
ple be large, even if the input-state map is surjective, 
resulting in large numerical errors for finite step sizes. 
Since the calculation of these higher derivatives is very 
expensive, we do not employ them in the calculation of 
the algorithmic step. The hypothesis that tracking glob- 
ally optimal paths in U(N) is typically more efficient 
than local optimization of the expectation value of the 
target observable is equivalent to the assumption that 
the second and higher order functional derivatives of the 
unitary propagator with respect to the control field are 
relatively small, but not negligible when attempting to 
traverse a large distance in U{N) in a single control field 
iteration. 

The computational expense of unitary matrix track- 
ing increases fairly steeply with system dimension. Since 
matrix inversion scales as N 2 , where N is the dimension 
of the matrix, the cost of inverting the G matrix scales 
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as TV 4 , where N is the Hilbert space dimension. By con- 
trast, global observable expectation value tracking, dis- 
cussed in the next section, avoids this overhead, but at 
the cost of being unable to specify precisely the unitary 
path followed during optimization. 

In order to test the hypothesis that the primary de- 
terminant of optimization efficiency is the unitary path 
length to the target W, we compared the optimization 
efficiencies of algorithms that follow a geodesic on the 
unitary group versus a faster path on the domain of ob- 
jective function values that corresponds to a longer path 
in U{N). These results are presented in section IVIII1 



IV. ORTHOGONAL OBSERVATION-ASSISTED 
QUANTUM CONTROL 

Unitary matrix tracking has the distinct advantage 
that it can directly follow an optimal path in the space 
of unitary propagators, assuming the input-state map 
is surjective and the linear formulation of the tracking 
equations above is a reasonable approximation. How- 
ever, it cannot be implemented experimentally without 
expensive tomography measurements, and carries a com- 
putational overhead that scales exponentially with sys- 
tem size. 

Given the initial state p(0) of the system, matrix 
elements of the unitary operator U(T) can be deter- 
mined based on knowledge of the final state p(T) = 
U(T)p(0)W{T). p{T) can be known only if a so-called 
tomographically complete set of observables has been 
measured sufficiently many times on identical copies of 
the system to approximate the expectation value of each 
observable. Assuming p(0) is nondegenerate, if such 
measurements are made at each step of the control op- 
timization, the unitary matrix tracking described above 
can be implemented (if p(0) is degenerate, the maximum 
number of U (T) elements that can be reconstructed will 
be diminished, as shown below). However, the cost of 
this procedure is very steep for large systems. The natu- 
ral question arises as to what sort of comparative benefit 
in optimization efficiency can be accrued from measure- 
ment of a limited number m of (orthogonal) operators, 
where m < iV 2 . 

Consider the case where n > m distinct observables, 
denoted 9i(T), 0„(T) or {6fc}, possibly linearly de- 
pendent and not necessarily orthogonal, are measured 
at each step. For simplicity represent each of the Her- 
mitian matrices as an iV 2 -dimensional vector with real 
coefficients. Then by Gram-Schmidt orthogonalization, 
it is always possible to construct an orthogonal ba- 
sis of m linearly independent iV 2 -dimensional vectors, 
&i(T), Q' m (T) that (spans this set) - any element 
of the set {Ok} can be expressed as a linear combina- 
tion of the basis operators in this set, i.e., for any k, 
8fe = J2iLi CikQ'i- In other words, the information ob- 
tained by measuring the expectation values of the set 
{Ok} is equivalent to that obtained by measuring the 



{9-}, since for each k, (8 fc ) = (E"=i c ifc©i)- orthogo- 
nal bases of Hermitian operators are the Pauli (2-d) and 
Gell-Mann (3-d) matrices. 

As above, we restrict ourselves here to coherent quan- 
tum dynamics, and additionally assume that a sufficient 
number of measurements of each observable have been 
made to accurately estimate its corresponding expec- 
tation value. Now consider the m (scalar functions of 
algorithmic time) {(Gfc(T, s))} of expectation values for 
each observable corresponding to a desired unitary track 
U(T,s). Again, the information about the states of the 
system p(T, s) or equivalently, U(T, s) contained in these 
measurements is equivalent to that contained in the m 
functions {(0J(T, s))}. Let us therefore represent this 
information in the form of the m-dimcnsional vector 
v(T, s), where 



Vi(T,s) 



5>*ei(r, a )>. 



During control optimization, we are interested in track- 
ing these paths v(T, s) in the vector space V that are 
consistent with the desired path Q(T,s) in U(N). 

The generalized differential equation (analogous to 
eqn ) that must be satisfied in order to simulta- 
neously track these paths is: 



dv(T, a) _ r T MT, a) gfM) dt = 
ds J 5e(s,t) ds 

f T / m^lflfv n \dQ(7>) 



(19) 



Based on eqn J2]), we have 



5e(s,t) ih 



-Tr 



^c fe4 e,(T),p(o) 



for the gradient of each of the observable expectation 
values Following the above derivation, we can con- 

vert this generalized differential equation into a vector 
integral equation: 




^ Cfci 6;(T),p(0) 



ds 



dQ(T,s) 
ds 



(20) 



Denoting the vector observable track of interest by w(s), 



i.e. 



w(s) = J2 Tl { P(°)Q f (T, s) c ^ 9 * Q( T ' s ) e » 



s 



and expanding d £ g*^ on the basis of orthogonal observ- 
ables, 



while 



de(s, t) 
ds 



6vi(T,a) 
Se(s,t) 



we have 



, T ( Sv(T,a) 
o V Se(s,t) 



or equivalently, 



E 

3=1 



(I 



Jvj(r,s) 

6s(s, t) 



oe{s, t) 



6e(s,t) 



dw(s) 



dw(s) 
ds 



Defining the correlation matrix in this case as 
■" T <5v,(7» 5vj(T, s) 



Tij (s) 



5s(s,t) Se(s,t) 



-dt, 



we obtain the following nonsingular algebraic differential 
equation for the algorithmic step in the control field: 



dE 



fs(t) 



dw 



-a(s) 



Se(s, t) 



(21) 



where f s (t) is again a free function and we have defined 
the vector function a(s) by analogy to a(s) above: 



a( S ) 



6v(T,s) 
Se(s, t) 



fs(t)dt. 



The advantage of orthogonal observable expectation 
value tracking, compared to unitary matrix tracking, is 
that the likelihood of the matrix T being ill-conditioned 
- even at abnormal extremal control fields s(t), where 
G is singular - diminishes rapidly with N 2 — m, where 
m is the number of orthogonal observable operators em- 
ployed. 

In the special case where only the observable of inter- 
est 8i is measured at each algorithmic step, this equa- 
tion reduces to: 



de(s, t) 
ds 



= f(s,t) + 



dP 
ds 



Jo a (s,t,T)f(s,t)dt 



7(s) 



a (s,t), 
(22) 

where P(s) is the desired track for (6i(T)}, a (s, t, T) = 
-^Tr (p(0) [U\T, 0)6 1 [/(T, 0), I7+(t, 0)fi(s)U(t, 0)] ) , 

and 7(s) = J [ao(s, t, T)] 2 dt. Here, it is of course 
not necessary to carry out any observable operator 
orthogonalization. 

Of course, measuring the expectation values (or gradi- 
ents) of two or more observable operators is more expen- 
sive than following the gradient of a single observable. 



However, note that the gradients ;; and — — - 



mSDl and ^MZ1> 

oe(tj 

of multiple observables are closely related since 
S(Qx) 



s = -^Tr{[C7t( T ) 0lf /( T ) )M ( t )] p(0 )} 



6(02) 
Se(t) 



= -jTr{[uHT)0 2 U(T)^(t)] p(0)}. 



As such, the information gathered through the estima- 
tion of the gradient of (8i(T)) can be used to "in- 
form" the estimation of (Oi(T)). In particular, al- 
though the norms of these two gradients differ, their 
time-dependencies - i.e., 5 fe7^yV~ fe7tfr^ are identical. 
Hence, only one high-dimensional gradient estimation 
needs to be carried out. 

The above algorithm can be applied to follow an 
arbitrary set of observable expectation value tracks 
{(0;(s))}. Here, we are interested in following the ob- 
servable tracks that correspond to the shortest path be- 
tween Uq and W on the domain of unitary propagators, 
namely the geodesic path U(s) = U(0) exp(i^4s) with 
A = — ilog(W*U(s)). As mentioned, the matrix W can 
be determined numerically if p(0) and arc known, for 
minimal computational cost. As shown by Hsieh et al. 
[3| , there exists a continuous submanifold of unitary ma- 
trices W that solve the observable maximization prob- 
lem; if we denote the Hilbert space dimension by N, the 
dimension of this submanifold ranges from N in the case 
that p(0) and are full rank nondegenerate matrices to 
N 2 — 2N + 2 in the case that p and are both pure 
state projectors (see section |VH| . 

The dimension of the subspace Mt of U(N) that is 
consistent with the observed track dv ^' s ' ) displays a 
complicated dependence on the eigenvalue spectra of 
p(0) and {0i}. We demonstrate this explicitly for the 
case of single observable tracking. In this case, 

M T = {V(s) | Hr(V(a)tp(0)V(a)e) = 

Tr(t/( S )tp(0)E/( S )e) = (0( S )» (23) 

where U(s) = exp (i log (W^Uq)s). As a first step, we 
must characterize the degenerate subset M(s) of uni- 
tary matrices that are compatible with a given ob- 
servable expectation value (0), as a function of the 
eigenvalue spectra of p(0) and 0. Let p(0) — Q^eQ 
and = ffi XR, where ei,£2,--- and Ai, A2, ... are 
the eigenvalues of p(0) and with associated unitary 
diagonalization transformations Q and R respectively. 
Then the observable expectation value corresponding to 
a given unitary propagator can be written 

J{U) = Tr(J7 t i?p(O)i? t ?7505 t ) 

= Tr:\{R^US) ] p{rfUS)Q\ = Tr(U^p(0)UQ) (24) 

where the isomorphism U — R^US also runs over U (N), 
and U = U(s). Therefore, without loss of generality, 
we can always assume that both p(0) and are in 
diagonal form, and determine M(s) in terms of U{s), 
instead of U(s). Denote by W(n) the product group 
lA(ni) x • ■ • x U(n r ), where U(ni) is the unitary group 
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acting on the n^-dimensional degenerate subspace corre- 
sponding to Xi, and define U(m) — U(mi) x ■ ■ ■U(m s ) 
in the same manner. Then any transformation U — > 
QUT\ where Q € U(n) and T G W(m), leaves J invari- 
ant: 

Tr(TC> t Q t ( 5(0)g?7T t e) = Tr([/ t Q t p(0)Q[/T t eT) 

= J(fi). 

This can be seen by observing that since p(0) is diag- 
onal, any unitary transformation p(0) — > Q'p(0)Q = 
p(0)Q'Q = p(0), if the unitary blocks of Q are aligned 
with the degeneracies of p(0). By the cyclic invariance 
of the trace, we also have 6 -» T^QT = QT^T = 8. 
Thus, the degenerate manifold can be written M(s) = 
U{n)U(s)U(m). Hence, the entire subspace of U(N) 
that is accessible to the system propagator during global 
observable tracking is 

M T = |J W(n)E7(a)W(m). 

0<s<l 

The manifold M(s) can be expressed as the quotient 
set 

M(s) = ^ (n)X ^ (m) . , 
U(m)nW(s)U(n)U(s) 

which can be seen as follows. Define Fjj{P,Q) ■ H — > 
PffQ, where H G W(iV) and (P, Q) G W(n) x W(m). 
Let stab(H) denote the stabilizer of H in W(n) x U(m), 
i.e. the set of matrix pairs (X,Y) G ZY(n) x t/(m) such 
that = XiJF = if. The stabilizer character- 

izes the set of points that are equivalent with H , hence 
the manifold M can be identified as the quotient set 
of W(n) x W(m) divided by stab(H). We can specify 
the stabilizer as follows. First, from Y = H^U^H, we 
see that H transforms U G U(n) into W(m). Hence 
y G W(m) n H^U(n)H. Conversely, for any V G 
W(m) r\H^U{n)H, the pair (Y+ffF, F) must be a mem- 
ber of stab(H). Hence, the stabilizer is isomorphic to 
U(m) nH^U(n)H. In the present case where H = U(s), 
we have 

stab(U(s)) = {(U{s)Y*U*(s),Y) : 

Y G W(m) n?7 t (s)W(n)t>(s)}. (25) 

Thus the dimension of the degenerate manifold M(s) is 

D (M(s)) = dimW(n) + dimW(m) - dim stab(U(s)). 

The dimension of this subspace cannot be specified 
in a simple form for arbitrary U(s) G U(N), since it is 
governed by the dimension of the stabilizer. We note 
a couple of special cases. If U(s) = R^U(s)S contains 
unitary subblocks that fall within the overlapping uni- 
tary subblocks in U(n) and £/(m),the dimension of the 
stabilizer is at least as large as that of the subblocks. If 
U (s) = BJU(s)S is a permutation matrix n or a product 



of a permutation matrix with a matrix of the aforemen- 
tioned type, the permutation matrix can act to rearrange 
the subblocks of U(m) so that they overlap with those 
of U(n) and thereby increase the dimension of the sta- 
bilizer. The dimension of the subspace (subgroup) of 
U{N) composed of such matrices U(s) can be shown to 
increase very rapidly with increasing degeneracies m and 
n in p(0) and O, respectively. 

The matrix U (s) can only rearrange or diminish the 
size of existing subgroups U(rii) within lA(n), but can- 
not create larger subgroups. Thus, we can establish a 
maximal dimension for stab(U(s)) for any given U(s) 
as that which maximizes the overlap between the re- 
spective subgroups U(rrij) and U(rii). This bound is 
achieved when the conjugation action U(s)'U(n)U(s) 
of U(s) is equivalent to the action n(s)tW(n)n of the 
permutation matrix n that rearranges the subblocks 
such that they display maximal overlap. Therefore, for 
fixed /o(0), O, the dimension of the manifold can range 
from the maximal value of nf 4- Y] j raj — N down to 

Ei=l n i + Ej=l m ) - El<Kr,l<j< a k % in the CaSe that 

the conjugation action of U(s) satisfies the above con- 
dition, where kij denotes the number of positions in the 
diagonal where the eigenvalues and ej appear simul- 
taneously after imposition of the permutation matrix n. 

Thus, we see that the size of the set of unitary matri- 
ces V producing the same observable expectation value 
(O) will change along the trajectory U(s), based on the 
extent to which the latter reorients the eigenvalues of 
p(Q) and 9 such that degeneracies coincide. In partic- 
ular, the volume of the subspace Mt of U{N) that is 
consistent with a track (0(s)) derived from a geodesic 
between Uq and W will display a strong dependence on 
the choice of Uq and W, and will change depending on 
the matrix W that is chosen from the degenerate sub- 
manifold of unitary matrices that solves the observable 
maximization problem. 

Note that the coefficient ao in the (single) observable 
expectation value tracking differential equation is in fact 
equal to the gradient on the domain e(t), equation 
Recall that the gradient flow (|4]) is defined by the differ- 
ential equation = ^jw-- Since the coefficients of a 
are scalars, we see that the algorithmic path for scalar 
tracking can be expanded on a basis whose dimension is 
identical to that of the gradient basis ([3|), as expected. 
We will analyze the dependence of the dimension of this 
basis on the eigenvalue spectra of p(0) and in section 

eh 

As a function of the algorithmic step s, the coeffi- 
cient /o T MM,r)/( s ,t)dt . n ^ ad . ugt gtep 

direction so that unitary matrices V(s) at each step are 
constrained within the subspace M(s). According to 
the above analysis, this dimension of this subspace will 
scale more steeply with increasing degeneracies in p(0) 
and O. The maximal dimension of M(s) ranges from 
N 2 — N for the problem where p and are both pure 
state projectors, to N for the case where p and are full 
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rank with completely nondegenerate eigenvalues. Even 
in the former case, this represents an advantage over the 
e(t)-gradient flow, which is free to explore the full N 2 - 
dimensional space of dynamical propagators in U{N). 

The above analysis assumes that the initial density 
matrix p(0) is known to arbitrary precision. This infor- 
mation is, of course, not required for e(i)-gradient based 
optimization, but it is readily acquired in the case that 
the initial state is at thermal equilibrium. For orthog- 
onal observable tracking, the cost of partial quantum 
state reconstruction of p(T) at each algorithmic step 
must be weighed against the increase in efficiency ob- 
tained by virtue of following a globally optimal path 
(section UX|) . 



apply error correction algorithms that exploit the curved 
geometry of the manifold. We therefore choose the error 
correction term to be a simple scalar multiple of the 
difference between the current values of the observable 
vector and its target value, i.e. (3 [w(s) — v(s)], such 
that the tracking differential equation becomes 



OS 



dw 



1 T 



f3 (w(s) - v(s)) + -^~ a(s) 



For the special case of single observable tracking, this 
reduces to 



V. ERROR CORRECTION AND FLUENCE 
MINIMIZATION 

1. Error correction 

In attempting to track paths on U(N), errors will in- 
evitably occur for two reasons. First, the algorithmic 
step on IA(N) will be a linear approximation to the true 
increment SU (T) due to discretization error; this error 
will increase as a function of the curvature of the inte- 
grated flow trajectory at algorithmic time s. Second, the 
D-MORPH integral equation is formulated in terms of 
only the first-order functional derivative 

^^Sp- for orthogonal observable and single observable 
tracking, respectively); the error incurred by neglecting 
higher order terms in the Taylor expansion will depend 
on the system Hamiltonian. 

In our numerical simulations, we apply error- 
correction methods to account for these deviations from 
the track of interest. (These methods can in principle 
also be implemented in an experimental setting.) For 
unitary matrix tracking, we correct for these inaccura- 
cies by following the (minimal-length) geodesic from the 
real point U (T, Sk) to the track point Q(sfc). This correc- 
tion can be implemented by incorporating the function 
■ \og(Q jt (sk)U(T, Sk)) into the matrix 



de(s, t) 
ds 



= f(s,t) + 



8e(t) ( m Se(t) 



differential equation for the algorithmic time step: 

-q- s = f(s, t) + (v(C(s)) + v(A) - a) G-y M (*)) 

In a more efficient approach, we combine error 
correction and the next gradient step in one it- 
eration [TB|. In this case, we define A(sk) = 
-^77=^ log {QHs k+1 )U(T, s k j) and use 

— = f(s > t)+(v(A)-a) G~ 1 v([i(t)). 

For orthogonal observable expectation value tracking, 
the vector space within which v(s) resides is not a Lie 
group, and consequently it is not as straightforward to 



/3(P(s) - (6(g)) + - a (s,t,T)f( S ,t)dt 

7(s) 



a (s,t). 
(27) 



2. Fluence minimization 



Clearly, the above analysis does not take into account 
the common physical constraint of penalties on the total 
field fluence. The effect of the fluence penalty in this sce- 
nario is then to decrease the degeneracy in the solutions 
to the above system of equations for d6 Q S f^ ■ This is ac- 
complished by choosing the free function /(s, t) in cither 
the unitary or orthogonal observable tracking differen- 
tial equations to be an explicit function of the electric 
field. It can be shown [ll[ that the choice: 

f(s,t) = -^e(s,t)W(t), 
As 

where W(t) is an arbitrary weight function and the As 
term controls numerical instabilities, will determine the 
d*^'*- 1 at each algorithmic time step s that minimizes 
fluence. 



VI. INTEGRATION OF QUANTUM 
OBSERVABLE EXPECTATION VALUE 
GRADIENT FLOWS 

We have seen (eqn ^ that in general, the projected 
path in U{N) that originates from following the lo- 
cal e(t)-gradient depends on the Hamiltonian of the 
system through the matrix G. Nonetheless, there is 
still a Hamiltonian independent component to the s(t)- 
gradient, which corresponds to the gradient on the do- 
main U(N). Thus, it is of interest to gain some un- 
derstanding of the behavior of this [/-gradient flow - 
whether it follows a direct path toward the target uni- 
tary propagator, or whether it biases the e(t)-gradient 
flow to follow indirect paths in U{N). Such an analysis 
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may shed light on the comparative optimization efficien- 
cies of the gradient compared to the tracking algorithms 
described in the previous sections. 

It can be shown that the U(T) and e(t) gradient flows 
evolve locally on subspaces of the same dimension, and 
that this dimension changes predictably as a function of 
the eigenvalue spectra of p(0) and 6. These gradient 
flows evolve on a subspace of the homogeneous space 
of U(N) whose dimension is given by the spectrum of 
the initial density matrix p(0), necessitating the use of a 
distinct coordinate basis to express the integrated gra- 
dient flow trajectories for different classes of p(0) that 
depend on the latter's number of nonzero and degener- 
ate eigenvalues pjj- Specifically let p(0) consist of r 



J 



subsets of degenerate eigenvalues pi, 
tiplicities n\,--- ,n r . Writing p(0) = 
have 



■ • ,p r , with mul- 
E"=iPiK)(*1> we 



<5$ 



i=l 



= a (t,T) 



N 

E 



(i|6(T)|i)(i|/x(t)|i)-(i| M (t)|j)(i|e(T)|i)] 

(28) 



Defining Sk 
written 



i=l ' 



k = 2 r 



1 this can be 



= E 



Se(t) 



k—1 i— Sfc + 1 

r 8k+i 

E 

fc=l Z— Sfc + 1 



s fc Sfc+i iV 

E+ E + E 



E+ E 

i=i i=s fc+ i+i 



{<i|e(r)li>0'K*M - (i|M*)li)(3'|e(r)|i)} (29) 
! /|e(T)|i)0Xi)|i) - mt)\j){j\e(T)\i)} (30) 



where the second equality follows from the fact that 
the terms corresponding to Ylj=^s h +i (i- e -j those aris- 
ing from the same degenerate eigenvalue of p(0))) are 
zero. This indicates that the dimension of the subspace 
of the space of skew-Hermitian matrices upon which the 
gradient flow evolves is 

r r 

D = N 2 - {N -nf n\ = n{2N - n) - ^ n\. 

This is the dimension of a compact polytope P which is 
the convex hull of the equilibria of the gradient vector 
field. The gradient flow involves on the interior of this 
polytope [181 ]. 

Both the e(t)-gradient flow (|4]) and observable track- 
ing (|2"2"]) can be expanded on this basis corresponding to 
j|^y . It can be shown (see below) that in the case of the 
gradient, the increased dimension of this basis set for 
increasing nondegeneracies in p(0), 9 generally results 
in increased unitary pathlengths; since the entire uni- 
tary group is free for exploration, the increased number 
of locally accessible directions results in the path mean- 
dering to more distant regions of the search space. By 
contrast, in the case of observable tracking, the increased 
number of locally accessible directions for greater non- 
degeneracies in p(0), are coupled with a decrease in 
the dimension of the globally accessible search space in 
U{N); the greater local freedom is used to follow the 
target unitary tracks of interest. 

In order to shed light on the origin of the aforemen- 
tioned behavior of the e(i)-gradient, we consider the 



global paths followed by [/-gradient flow of $. In a use- 
ful analogy, the control optimization process can itself 
be treated as a dynamical system. The critical manifolds 
of the objective function then correspond to equilibria of 
the dynamical system, and the gradient flow trajectories 
to its phase trajectories. Within this analogy, the gra- 
dient flow of <I> on U{N) can be shown to represent the 
equations of motion of an integrable dynamical system. 
The expression ([5]) for the gradient flow of $ above is 
cubic in U. However, through the change of variables 
U(s, T) -> p(s, T) = U(s, T)p(0, 0)W(s, T) we can reex- 
press it as a quadratic function: 

p(s,T) = 

- U(s,T)p(Q t 0)U*(s,T) - U(s,T)p(0,0)U i (s,T) 
- p 2 (s, T)G - 2p{s, T)@p(s, T) + Gp 2 (s, T) 

= [p(s,T),[p(s,T),e]} (31) 

where s denotes the algorithmic time variable of the gra- 
dient flow in U(N) and the dot denotes the s-derivative. 
This quadratic expression for thegradient flow is in so- 
called double bracket form 

USEE El- The set of all 

U(s, T)p(0, 0)C/ t (s, T) is a homogeneous space M(p) for 
the Lie group U(N), namely the space of all Hermi- 
tian matrices with eigenvalues fixed to those of p(0,0). 
Maximizing $(J7) over IA(N) is equivalent to minimiz- 
ing the least squares distance ||0 — UpU^\\ 2 of to 
UpW E M{p): 

= ||0f -2$([/) + |H| 2 . 
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Here, we provide an explicit formula pertaining to 
the analytical solution for the above gradient flow for 
what is perhaps the most common objective in quan- 
tum optimal control theory and experiments, namely the 
maximization of an arbitrary observable starting from a 
pure state. In particular, this includes the special case 
of maximizing the transition probability Pif between 
given initial and final pure states \i) and |/). When- 
ever p(t — 0) is a pure state, n = r = ri\ = 1, and 
D = 2N — 2, and the convex hull of the critical points of 
the vector field is a (N — l)-dimcnsional simplex. The 
gradient flow evolves in the interior of this simplex. The 
analytical solution for the fully general case of a mixed 
state p(0) and nondegenerate is more complicated and 
presented in another work of the authors. 

Since the objective function is symmetric with respect 
to p(0) and 0, this formulation applies if either p(0) or 
O is a pure state projection operator, i.e., if at least one 
of them can be diagonalized by an appropriate change 
of basis to matrices that have only one nonzero diagonal 
element, corresponding to or |/)(/|, respectively. 

The other operator can have an arbitrary spectrum. The 
same integrated gradient flow thus applies to the prob- 
lem of maximizing the transition probability between 
any generic mixed initial state to any pure state. 



Under these conditions, we can execute a change 
of variables such that the double bracket flow, which 
evolves on the ^m(m + l)-dimensional vector space of 
Hermitian matrices p(0) is mapped to a flow on the m- 
dimensional Hilbert space. Letting \ip(s)) = U(s)\i), the 
double bracket flow can be written: 

W*)> = U(»)\i) 

= [QU(s)\i) -U(s)\i)(i\uHs)QU(s)]\i) 

= [e-<v>00|ew*)>J] M«)>. 

If we define x(s) = (\ci(s)\ 2 , ■ ■ ■ , \cn(s)\ 2 ), where 
Ci(s), • • • , cjv(s) are the coordinates of \ip(s)) under the 
basis that diagonalizes 0; it can be verified that the in- 
tegrated gradient flow can be written: 



o2se 



x(s) = 



M0)| 2 



EtiK(0)| 2 e 
| Cl (0)|V-- ,e 2 ^«M0)| 



2 e 2sXt 



a 2s\ 1 



Eli\cM\ 2 e 



2sXi 



(32) 



(33) 



where Ai, ■ • ■ , Ajv denote the eigenvalues of 0. 

In the case that has only one nonzero eigenvalue 
(pure state), this becomes: 



x(s) = 



2sX 



|ci(0)| s 
E^|ci(0)| 2 +e 

|ci(0)| 2 
l + (e2«Aj -l)| Cj (0)| 2 ' 



^2s\j 



Mo) | 5 



Mo)| 2 '' 



E£>(0)| a + e^| Cj (0)|i 



c,(0)| s 



l + (e 2 ^ -l)| Cj (0)| 2 ' 
I 



In contrast to quantum time evolution of the state 
vector (which resides on the Hilbert sphere S 1 ™ -1 , the 
matrix e 2 translates the vector x, which resides on the 
(m — l)-dimensional simplex, through algorithmic time. 
Since coherent quantum time evolution cannot change 
the eigenvalues of the density matrix p, the optimization 
of controls simply reorders these eigenvalues. Hence the 
gradient flow is said to be an " isospectral" flow. We are 
primarily interested in the gradient flow on the domain 
of unitary propagators of the quantum dynamics, U(T). 
In general, there exists a one-to-many map between p(T) 
and U(T). This attests to the existence of a multiplic- 
ity of paths through the search space of the quantum 
optimal control problem that will maximize the objec- 
tive functional in equivalent dynamical time. Although 
the analytical solution we have presented is framed on 
the homogeneous of U(N) rather than on the Lie group 
itself, we will see in the next section that generic proper- 
ties of the behavior of the flow trajectories, in particular 
their phase behavior with respect to local critical points, 
are identical on these domains. 

The above gradient flow for only applies when 



has no degeneracy in its eigenvalues, including zero 
eigenvalues. If the maximum eigenvalue of is degen- 
erate with multiplicity fc, such that c\ = c% = ■ ■ ■ = Cfc, 
and Ck > Cj, j = k + 1, • • • , n, then the dynamics con- 
verges to the point |(1, • • • , 1, 0, ■ ■ • ,0). 

A remarkable feature of the gradient flow for objective 
function $ is that it is a Hamiltonian flow 0, [2l|, , 
for general p(0) and 0. The eigenvalues can be viewed as 
the analog of the momenta in the corresponding Hamil- 
tonian system. The "isospectral" character of the flow 
indicates these momenta are conserved. An alternative 
proof of the integrability of the flow for $ is based on 
demonstrating that in N dimensions, the flow has N in- 
tegrals of the motion that are in involution, which is the 
classical definition of complete integrability for a Hamil- 
tonian system. From the point of view of the modern 
theory of integrable systems, the double bracket flow 
can be shown to represent a type of Lax pair, a general 
form that can be adopted by all completely integrable 
Hamiltonian systems [23j. 
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VII. PHASE BEHAVIOR OF QUANTUM 
OBSERVABLE MAXIMIZATION GRADIENT 
FLOWS 

The integrated flow trajectories provided above for the 
gradient of $ on the domain of unitary propagators can 
be used to provide insight into global behavior of the 
e(i)-gradient flow. Just as the global trajectory of the 
unitary geodesic flow influences, but does not directly 
determine the unitary path followed by global observ- 
able tracking, the integrated [/(T)-flow trajectories in- 
fluence but do not completely determine the behavior 
of the e(i)-gradient flow. A useful metric for assessing 
the phase behavior of the i7(T)-gradient flows is the dis- 
tance of the search trajectory from the critical manifolds 
of the objective as a function of algorithmic time. 

The distance of the search trajectory to the global 
optimum of the objective can be expressed: 



\\ X (s)-e l 4 2 = Ms)\\ 2 -2 



(e^x(0),e„) 
(e 2se x(0)) 



The equilibrium points of this flow are the critical points 
of the objective function. Hsich et al. [|[ showed that 
the critical manifolds of <& satisfy the condition 



d,/ 



iTr(A [e,£/V(0)£7]) = 0. 



The critical manifolds are then given by matrices of the 
form 

Ui = QPiR\ 

where Pi, I = 1, • • • , N\ is an AT- fold permutation 
matrix whose nonzero entries are complex numbers 
exp(i^i), exp(i^jv) of unit modulus, and p = Q^eQ 
and 9 = R^XR. The critical manifolds have a similar 
topology to the manifolds M(s) shown in section ITVl to 
map to a given observable expectation value. For exam- 
ple, in the case that p and are fully nondegenerate, 
they are ./V-torii TJ . In this case, the number of crit- 
ical manifolds scales factorially with the Hilbert space 
dimension N. In the case that either p or O has only 
one nonzero eigenvalue, while the other is full rank and 
nondegenerate, the number of critical manifolds scales 
linearly with N. (In the present case, assuming is 
full-rank, the number of these equilibrium points scales 
exponentially as 2 m with the Hilbert space dimension.) 
The optimal solution to the search problem corresponds 
to the basis vector where has its maximal eigenvalue. 
In the case that the observable operator has only one 
nonzero eigenvalue i, there are only two critical points, 
corresponding to ±e*. 

The time derivative of the distance of the search tra- 
jectory to a critical manifold (framed on the homoge- 
neous space) is 



dt 



11*00- 



(o)EL(^-^r Aj *<(o) 



(34) 

Solving for the zeroes of this time derivative reveals that 
the distance between the current point on the search 
trajectory and the solution can alternately increase and 
decrease with time. The distance of the gradient flow 
trajectory from the suboptimal critical manifolds can 
alternately increase and decrease arbitrarily depending 
on the spectral structure of 0. Thus the density matrix 
does not always resemble the target observable operator 
to a progressively greater extent during the algorithmic 
time evolution. 

Qualitatively, the most obvious and important feature 
of these trajectories is that the closer the initial condi- 
tion is to a suboptimal critical manifold, the greater the 
extent to which the gradient flow follows the boundary of 
the simplex during its early time evolution. Away from 
these initial conditions, the behavior of the gradient flow 
trajectory is considerably more sensitive to the spectral 
structure of the observable operator than it is to the 
initial state f/'o ■ Indeed, it may be shown [24| that for op- 
erators with two eigenvalues arbitrarily close to each 
other, the time required for convergence to the global 
optimum increases without bound, whereas this is not 
the case for initial states with a arbitrarily close to 0. 

Of course, the search trajectory slows down in the 
vicinity of the critical manifolds. Note that at the criti- 
cal manifolds themselves, the observable tracking equa- 
tions also encounter singularities, since ao = 0. In the 
vicinity of the critical manifolds (the attracting regions) 
tracking may be less accurate due to small values of 'y(s) 
in the denominator of eqn (|22[) . corresponding to func- 
tions ao of small modulus. 

As the nondegeneracies in p(0) and increase, the 
number of attractors of the search trajectory increases, 
indicating an increase in the unitary path length. There- 
fore, as mentioned above, although the trajectory fol- 
lowed by the gradient locally accesses the same number 
of directions as observable tracking, its global path ap- 
pears biased towards being longer, assuming the Hamil- 
tonian is fixed, especially for highly nondegenerate p(0), 
0; the local steps in global geodesic observable tracking 
access the same number of local directions but orient 
them towards unitary matrices along a shorter path. 



VIII. NUMERICAL IMPLEMENTATION 

Simulations comparing the efficiency of unitary or or- 
thogonal observable tracking with gradient-based opti- 
mal control algorithms will be presented in a separate 
paper. However, we provide here a brief summary of 
numerical methods that can be used to implement the 
various tracking algorithms described above. 
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For multiple observation-assisted tracking, a set of 
n observable operators ©i,--< ,6 n were either chosen 
randomly or based on eigenvalue degeneracies. This 
(possibly linearly dependent) set was orthogonalized 
by Gram-Schmidt orthogonalization, resulting in a m 
(where m < N) dimensional basis set of observable op- 
erators. The m-dimensional vector v(s) was constructed 
by tracing each of these observable operators with the 
density matrix. 

Numerical solution of the D-MORPH differential 
equations (fl~5l I2T1 was carried out as follows. The 
electric field e(s, t) was stored as a p x q matrix, where 
p and q are the number of discretization steps of the 
algorithmic time parameter s and the dynamical time 
t, respectively. For each algorithmic step Sk, the field 
was represented as a g-vector for the purpose of com- 
putations. Starting from an initial guess e(so,t) for the 
control field, the Hamiltonian was integrated over the in- 
terval [0, T ] by propagating the Schrodinger equation 
over each time step tk — > ifc+i, producing the local 
propagator U(tj+i,tj) = cxp [—iH{si,tj)T/{q— 1) ]. 
For this purpose, the propagation toolkit was used. Lo- 
cal propagators were precalculated via diagonalization 
of the Hamiltonian matrix (at a cost of N 3 ), exponen- 
tiation of the diagonal elements, and left/right multipli- 
cation of the resulting matrix by the matrix of eigenvec- 
tors/transpose of the matrix of eigenvectors. This ap- 
proach is generally faster than computing the matrix ex- 
ponential directly. Alternatively, a fourth-order Rungc- 
Kutta integrator, can be employed for the propagation, 
and is often used for density matrix propagation. 

The time propagators U(tj,0) = 
U(tj, ij-i), • • • , U(ti, to) computed in step 1 were 
then used to calculate the time-evolved dipole operators 
fi(tj) — W(tj,0)fiU(tj,0), which can be represented as 
a g-dimensional vector of N x N Hermitian matrices. 
The iV 2 x TV 2 matrix G(sfc) and iV 2 vector a(sk) (alter- 
natively, the m x m matrix T(sk) and m vector a(sfc))) 
were then computed by time integration of the dipole 
functions (and appropriate choice of function f(s,t) 
)described above. For tracking of unitary gradient 
flows, the next point Q(sk) on the target unitary track, 
necessary for the implementation of error correction, 
was calculated in one of two ways: i) numerically 
through Q(sk) = Q(sfc-i) * exp(— iA(sfc)ds), computed 
using a matrix exponential routine (because of the 
greater importance of speed vs accuracy for this step) 
or ii) analytically, through the integrated flow equation 
above. If error correction was employed, the matrix 
C(s fe ) = - Sk _ l Sk l log(t/ t (sfc)Q(sfc)) for unitary error 
correction was calculated by diagonalization of the 
unitary matrix followed by calculation of the scalar log- 
arithms of the diagonal elements, where the logarithms 
are restricted to lie of the principal axis. 

Next, the control field e(sk,t) was updated to 
e(sk+i,t). This step required the inversion of the 
TV 2 x iV 2 matrix G(sfc) or m x m matrix r(s/ c ), 
which was carried out using LU decomposition. 



The quantities G 1 (sfe), A(sfc), a(sk) (alternatively, 
r _1 (sfe), ^p-,a(sfe), for scalar tracking) and C(sk) were 

used to compute the (^-dimensional vector 8g g*'^ ■ One 
of two approaches was used to update the field: (i) 
a simple linear propagation scheme, i.e. e(s / r c+1 ,<) — 

e(sfc,t) -I- ds - gj 2 ^ , or (ii) a fourth-order Runge-Kutta 
integrator. Because the accuracy of tracking depended 
largely on the accuracy of this s-propagation step, and 
because only one s-propagation was carried out for each 
set of q time propagations, a more accurate (but expen- 
sive) fifth-order Runge-Kutta integrator was often used 
in this step. The updated control field s(sk+i, t) was 
then again used to the propagate the Schrodinger equa- 
tion. 

In the above numerical tracking scheme, the most 
computationally intensive step is the propagation of the 
Schrodinger equation. Although calculation of the ma- 
trix elements of G for unitary tracking (respectively, 
r for orthogonal observable tracking) and inversion of 
this matrix scale unfavorably with the dimension of the 
quantum system, this scaling is polynomial. Since re- 
maining on the target t/-track can play an important 
role in the global convergence of the algorithm, es- 
pecially where the local gradient follows a circuitous 
route, the additional expense incurred in accurate s- 
propagation may be well-warranted. 

IX. IMPLICATIONS FOR OPTIMAL 
CONTROL EXPERIMENTS (OCE) 

The majority of optimal control experiments (OCE) 
have been implemented with adaptive or genetic search 
algorithms, which only require measurement of the ex- 
pectation value of the observable for a given control field. 
Recently, gradient-following algorithms have been imple- 
mented in OCE studies, based on the observation that 
the control landscape is devoid of suboptimal traps [H . 
In computer simulations, it is generally observed that 
gradient-based algorithms converge more efficiently than 
genetic algorithms. However, more sophisticated local 
search algorithms, such as the Krotov or iterative algo- 
rithms often used in OCT, are difficult if not impossible 
to implement experimentally due to the extreme diffi- 
culty in measuring second derivatives of the expectation 
value in the presence of noise and decoherence. There- 
fore, global search algorithms which can be implemented 
on the basis of gradient information alone, such as or- 
thogonal observable tracking, are particularly attractive 
candidates for improving the search efficiency of OCE. 

Applying observable tracking experimentally requires 
a fairly simple extension of gradient methods that have 
already been applied. The gradient is determined 
through repeated measurements on identically prepared 
systems, to account for the impact of noise. Instead of 
following the path of steepest descent, the laser field is 
updated in a direction that in the linear approximation 
would produce the next observable track value 0(s+ds). 
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The assumption is that because this observable step is 
consistent with a short step on the domain of unitary 
propagators, the associated error in reaching this expec- 
tation value will be smaller than that associated with 
trying to move over the same fraction of the gradient 
flow trajectory. Error correction can be implemented 
using a method identical to that described above for nu- 
merical simulations. 

In the simplest incarnation of orthogonal observation- 
assisted quantum control, the full density matrix p(0) 
is estimated at the beginning of the control optimiza- 
tion. This need be done only once and hence adds only 
a fixed overhead to the experimental effort that does not 
add substantially to the scaling of algorithmic cost with 
system dimension. At subsequent steps during the op- 
timization, the experimenter can (possibly adaptively) 
decide how many orthogonal observations to make in or- 
der to estimate the final density matrix p(T), with the 
goal of keeping the unitary path as close as possible to 
the desired geodesic. The number of distinct observables 
that must be measured at each step (and hence roughly 
the total number of measurements) scales linearly with 
the number of estimated parameters of p(T). 

In order to properly compare the efficiencies of experi- 
mental gradient-following and global observable tracking 
methods, it is necessary to consider the expense associ- 
ated with reconstruction of the initial and final density 
matrices p(0) and p(T) in the latter case. A variety 
of different quantum statistical inference methods have 
been developed over the past several years for the esti- 
mation of density matrices on the basis of quantum ob- 
servations. Like the measurement of the gradient, these 
methods are based on multiple observations on identi- 
cally prepared copies of the system. If we consider M 
measurements on identically prepared copies, each mea- 
surement is described by a positive operator- valued mea- 
sure (POVM). Of interest to us here is the scaling of the 
number of measurements required to identify the den- 
sity matrix within a given precision, with respect to the 
Hilbert space dimension. In most reconstruction tech- 
niques, such as quantum tomography (25j , a matrix ele- 
ment of the quantum state is obtained by averaging its 
pattern function over data. In the averaging procedure, 
the matrix elements arc allowed to fluctuate statistically 
through negative values, resulting in large statistical er- 
rors. Recently, the method of maximal likelihood esti- 
mation (MLE) of quantum states has received increasing 
attention due to its greater accuracy 3 . Denoting by Fj 
the POVM corresponding to the i-th observation, the 
likelihood functional 



3 Other methods for state reconstruction, such as the maximum 
entropy method or Bayesian quantum state identification l2tl . 
can also be employed. 



N 

L(p) =l[Tr(p(0)Fi) 

i=l 

describes the probability of obtaining the set of observed 
outcomes for a given density matrix p. This likelihood 
functional is maximized over the set of density matrices. 
An effective parameterization of p is p(Q) — T^T, which 
guarantees positivity and Hermiticity, and the condition 
of unit trace is imposed via a Lagrange multiplier A, to 
give 

N 

L(f) = ^ In Tr (ftfFj) - ATr(f t f). 

i=l 

Standard numerical techniques, such as Newton- 
Raphson or downhill simplex algorithms, are used to 
search for the maximum over the N 2 parameters of the 
matrix T. The statistical uncertainty in density matrix 
estimates obtained via MLE can quantified by consider- 
ing the likelihood to function to represent a probability 
distribution on the space of density matrix elements - 
or, in the current parameterization, the N 2 parameters 
constituting the matrix T, denoted here by the vector 
t. In the limit of many measurements, this distribution 
approaches a Gaussian. The Fisher information matrix 
/ = Q t Q t , , which is the variance of the score function 
that is set to zero to obtain the MLE estimates, can 
then be used to quantify the uncertainties in the param- 
eters. Note that the constraint Tr(TTT') = 1 implies that 
the optimization trajectory maintains orthogonality to 

u = aTr g^ T ^ ■ Under these conditions it can be shown 
that the covariance matrix for t is given by 

v = i^ 1 - r l uu T r l ju T r l u. 

As such, the associated uncertainties on the density ma- 
trix elements as a function of the number of measure- 
ments can be determined for a given system based on 
computer simulations. 

Banaszek et al. [2?} applied MLE to both discrete 
and continuous quantum state estimation, comparing to 
state tomography. In the case of continuous variable 
states, the density matrix was, of course, truncated. For 
identical systems, MLE required orders of magnitude 
fewer measurements TV to reconstruct the state with the 
same accuracy. For example, only 50, 000 homodync 
data (compared to 10 s for tomography) were required to 
reconstruct the matrix for a single-mode radiation field. 
The number of required measurements was not highly 
sensitive to the truncation dimension, since adaptive 
techniques can be used to improve efficiency in higher 
dimensions. Only 500 measurements were required to 
reconstruct the density matrix of a discrete quantum 
system, a pair of spin-1/2 particles in the singlet state. 
Given that accurate estimation of the gradient requires 
a similar number of measurements if MLE is used 
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for state reconstruction, the additional algorithmic over- 
head for observable tracking is not limiting. Quantita- 
tive calculations of this overhead will be reported in a a 
forthcoming numerical study. 

X. DISCUSSION 

We have presented several global algorithms for the 
optimization of quantum observables, based on following 
globally optimal paths in the unitary group of dynami- 
cal propagators. The most versatile of these algorithms, 
orthogonal observable expectation value tracking, aims 
to simultaneously track a set of observable expectation 
value paths consistent with the unitary geodesic path to 
the target propagator. The performance of the latter 
has been compared theoretically to that of local gradi- 
ent algorithms. A follow-up paper will report numerical 
simulations comparing the efficiencies of the algorithms 
described herein, across various families of Hamiltoni- 
ans. 

Although the e(i)-gradicnt flow is always the locally 
optimal path, its projected path in U{N) is generally 
much longer than those that can be tracked by global 
algorithms. The latter often require fewer iterations for 
convergence. Of course, in order to assess the utility of 
[/-flow tracking as a practical alternative to local OCT 
algorithms, it is necessary to consider the computational 
overhead incurred in solving the system of tracking dif- 
ferential equations at each algorithmic step, which is on 
the order of A 4 -times more costly than solving scalar 
observable tracking equations. 

For the systems studied, the geodesic track in U(N) 
can typically be followed faithfully by matrix tracking 
algorithms, assuming the system is controllable on the 
entire unitary group. However, for many Hamiltonians, 
unitary matrix tracking can routinely encounter regions 
of the landscape where the G matrix is ill-conditioned. 
In contrast, tracking a vector v(s) of m orthogonal ob- 
servable expectation values corresponding to this uni- 
tary matrix track encounters such singularities more 
rarely, provided m < N. 

However, the number of observable expectation values 
tracked is not the only factor that affects the mean path- 
length in U (N) . The relationship between the basis set 
of operators (spanning the space of measured observ- 
ables) and the eigenvalue spectrum of p(0) affects the 
dimension and volume of the subspace of U(N) that is 
accessible to the search trajectory. Indeed, the compar- 
ative advantage of employing global observable tracking 
algorithms versus local gradient algorithms was found 
to depend on the spectra of the initial density matrix 
j0(0), as predicted based on a geometric analysis of the 
map between unitary propagators and associated ob- 



servable expectation values. In particular, measurement 
of the expectation value of a quantum observable pro- 
vides more information about the dynamical propagator 
U(T) when p(0) has fewer degenerate eigenvalues, since 
degeneracies produce symmetries that result in invari- 
ant subspaces over which the unitary propagator can 
vary without altering the observable. Thus, for an iden- 
tical number of measurements, the global path in U (N) 
can be tracked with greater precision for systems with 
nondegenerate p(0). 

Clearly, a very important question underlying the ef- 
ficiency of global OCT algorithms based on paths in the 
unitary group is whether the nonsingularity of G and 
the assumption of small higher-order functional deriva- 
tives remains valid for more general Hamiltonians be- 
yond those considered here. In principle, paths in hi (N) 
that are longer than the geodesic might be significantly 
easier to track if these assumptions were to break down, 
for particular systems. For systems where G is almost 
always close to singular, global tracking algorithms on 
U(N) may not be viable, even in the presence of er- 
ror correction. In these cases one would expect global 
observable tracking to be preferable to matrix tracking 
algorithms, since the system can "choose" which of the 
infinite number of degenerate paths in hi{N) it follows. 

This study, and the associated forthcoming numerical 
simulations, sets the stage for experimental testing of 
its prediction that global observable tracking algorithms 
will display advantages compared to the gradient. As 
discussed above, an important question for future study 
is how to implement global observable control algorithms 
experimentally, and whether nonideal conditions (noise, 
decoherence) in the laboratory will obscure some of its 
predicted advantages. In particular, the effective degen- 
eracies of p(Q) and O - e.g., whether the populations of 
the various levels in a mixed state are sufficiently high 
to permit accurate determination of the associated uni- 
tary propagators in the presence of noise - become very 
important. 

Besides the perceived advantage of global observable 
control algorithms, they may offer a means of assessing 
the search effort and search complexity inherent in quan- 
tum control problems in a universal manner, since they 
are more system independent than local gradient-based 
algorithms. As shown in a separate numerical study 
comparing these algorithms, over several families of re- 
lated Hamiltonians, the variance of the convergence time 
of the e(i)-gradicnt flow is significantly greater than that 
of global observable or unitary matrix tracking. Given 
that its convergence is also faster, unitary tracking may 
offer an approach to setting upper bounds on the scal- 
ing of the time required for quantum observable control 
optimizations, as a function of system size. 
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